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Abstract. In this paper we consider finite-dimensional constrained Hamiltonian 
systems of polynomial type. In order to compute the complete set of constraints and 
separate them into the first and second classes we apply the modern algorithmic 
methods of commutative algebra based on the use of Grobner bases. As it is shown, 
this makes the classical Dirac method fully algorithmic. The underlying algorithm 
implemented in Maple is presented and some illustrative examples are given. 



1 Introduction 

The generalized Hamiltonian formalism invented by Dirac Q for constrained 
systems has become a classical tool for investigation of gauge theories in 
physics |2[^,H , and a platform for numerical analysis of constrained mechani- 
cal systems H . Finite-dimensional constrained Hamiltonian systems are part 
of differential algebraic equations whose numerical analysis is of great re- 
search interest over last decade |^] because of importance for many applied 
areas, for instance, multi-body mechanics and molecular dynamics. 

In physics, the constrained systems are mainly of interest for purposes 
of quantization of gauge theories which play a fundamental role in mod- 
ern quantum field theory and elementary particle physics. Dirac devised his 
methods to study constrained Hamiltonian systems just for those quantiza- 
tion purposes. Having this in mind, he classified the constraints in the first 
and second classes. A first class constrained physical system possesses gauge 
invariance and its quantization requires gauge fixing whereas a second class 
constrained system does not need this. The effect of the second class con- 
straints may be reduced to a modification of a naive measure in the path 
integral. The presence of gauge degrees of freedom (first class constraints) 
indicates that the general solution of the system depends on arbitrary func- 
tions. Hence, the system is underdetermined. To eliminate unphysical gauge 
degrees of freedom one usually imposes gauge fixing conditions whereas for 
elimination of other unphysical degrees of freedom occurring because of the 
second class constraints, one can use the Dirac brackets [H§,0|. In some spe- 
cial cases one can explicitly eliminate the unphysical degrees of freedom |^ . 

* This work was supported in part by Russian Foundation for Basic Research, 
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Unlike physics, where constrained systems are singular, as they contain in- 
ternal constraints, mechanical systems are usually regular with externally 
imposed constraints Such a system is equivalent to a singular one whose 
Lagrangian is that of the regular system enlarged with a linear combination 
of the externally imposed constraints whose coefficients (multipliers) are to 
be treated as extra dynamical variables. The latter system may reveal extra 
constraints for the former system providing the consistency of its dynamics. 

Therefore, to investigate a constraint Hamiltonian system one has to de- 
tect all the constraints involved, and separate them, for physical models, into 
first and second classes. In his theory |^ Dirac gave the receipt for compu- 
tation of constraints which is widely known as Dirac algorithm, and it has 
been implemented in computer algebra software pO| . However, the Dirac ap- 
proach, as a method for computation of constraints, is not yet an algorithm. 
Even computation of the primary constraints, given a singular Lagrangian, is 
not generally algorithmic. Moreover, in generation of the secondary, tertiary, 
etc., constraints by the Dirac method one must verify if a certain function of 
the phase space variables vanishes on the constraint manifold. Generally, the 
latter problem is algorithmically unsolvable. Similarly, there are no general 
algorithmic schemes for separation of constraints into the first and second 
classes. In physical literature one can find quite a number of particular meth- 
ods developed for the constraint separation (see, for example, ||ll|,|l2j). But 
all of them have non-algorithmic defects. Thereby, being successfully applied 
to one constrained system, those methods may be failed for another system 
even of a similar type. 

In practice, many constrained physical and mechanical problems are de- 
scribed by polynomial Lagrangians that lead to polynomial Hamiltonians. 
In this case, as we show in the present paper, one can apply Grobner bases 
which nowadays have become the most universal algorithmic tool in commu- 
tative algebra Q and algebraic geometry The combination of the 
Dirac method with the Grobner bases technique makes the former fully al- 
gorithmic and, thereby, allows to compute the complete set of constraints. 
Moreover, the constraint separation is also done algorithmically. We show 
this and present the underlying algorithm which we call algorithm Dirac- 
Grobner. This algorithm has been implemented in Maple V Release 5, and 
we illustrate it by examples both from physics and mechanics. 



2 Dirac Method 

In this section we shortly describe the computational aspects of the Dirac 
approach to constrained finite-dimensional Hamiltonian systems . 
Let us start with a Lagrangian L{q,q) = L(qi, qj) (1 < i,j < n) as a function 
of the generalized coordinates qi and velocities q^^. If the Hessian d'^L/dqidqj 

^ We consider only autonomous systems, and there is no loss of generality since 
time t may be treated as an additional variable. 
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has the fuU rank r = n, then the system is regular and it has no internally 
hidden constraints. Otherwise, ii r < n, the Euler-Lagrange equations 

= (l<.<n) (1) 

with 

9L 

are singular or degenerate, as not all differential equations (|l|) are of the 
second order. There are just n — r such independent lower order equations. 
By the Legendre transformation |^ 

Hc{p,q) ^ PiQi - L, (3) 

we obtain the canonical Hamiltonian with momenta pi defined in (|^) . In the 
degenerate case there are primary constraints denoted by (j)a, which form the 
primary constraint manifold denoted by Sq 

So ■■ MP,q)^0 {l<a<n~r), (4) 

Thus, the dynamics of the system is determined only on the constraint man- 
ifold (^. To take this fact into account, Dirac defined the total Hamiltonian 

Ht= Hc+Ualpa (5) 

with multipliers Ua as arbitrary (non-specified) functions of the coordinates 
and momenta. The corresponding Hamiltonian equations determine the sys- 
tem dynamics together with the primary constraints 

qi^{Ht,qi\, Pi = {Ht,Pi}, cj)a{p,q)^Q {I < i < n, 1 < a < n - r), (6) 

where the Poisson brackets are defined for any two functions /, g of the dy- 
namical variables p and q as follows 

r . . df dg dg dp 
opi oqi opi oqi 

In order to be consistent with the system dynamics, the primary constraints 
must satisfy the conditions 

^a^{Ht,M = (l<a<n-r), (8) 

where stands for the equality, called a week equality, on the primary con- 
straint manifold (^. The Poisson bracket in (||) must be a linear combination 
of the constraint functions |^ . Given a constraint function 0q , the consistency 

^ In this paper summation over repeated indices is we always assumed. 
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condition (^), unless it is satisfied identically, may lead either to a contradic- 
tion or to a new constraint. The former case signals that the given Hamilto- 
nian system is inconsistent. In the latter case, if the new constraint does not 
involve any of multipliers Ua, it must be added to the constraint set, and, 
hence, the constraint manifold must involve this new constraint. Otherwise, 
the consistency condition is considered as defining the multipliers, and the 
constraint set is not enlarged with it. 

The iteration of this consistency check ends up with the complete set of 
constraints such that for every constraint in the set condition (H) is satisfied. 
This is the Dirac method of the constraint computation. As shown in |p^ , 
the method is nothing else than completion of the initial Hamiltonian system 
to involution, and the constraints generated are just the integrability condi- 
tions. For general systems of PDEs, the completion process is done by 
sequential prolongations and projections. For Hamiltonian systems, the time 
derivative of a constraint is its prolongation whereas projection of the pro- 
longed constraint is realized in (^) by computing the Poisson bracket on the 
constraint manifold. 

Let now S be the constraint manifold for the complete set of constraints 

U : 0„(p,g) = O (l<a<fc). (9) 
If a constraint function 0^ satisfies the condition 

{(t>a{p,q),MP.<l)}^0 (l</3<fc), (10) 

it is of the first class. Otherwise, the constraint function is of the second class. 
The number of the second class constrains is equal to rank of the following 
(fc X k) Poisson bracket matrix, whose elements must be evaluated on the 
constraint manifold 

= (11) 

Note that matrix M has even rank because of its skew-symmetry. 

If a Lagrangian system LQ{q,q) is regular with externally imposed holo- 
nomic constraints ipaio) — 0, the system is equivalent Q to the singular one 
with Lagrangian L ^ Lq + Xa4>a and extra generalized coordinates Aq. Fur- 
thermore, the Dirac method can be applied for finding the other constraints 
inherent in the initial regular system and, hence, not involving the extra 
dynamical variables. 

Therefore, the problem of constraint computation and separation is re- 
duced to manipulation with functions of the coordinates and momenta on the 
constraint manifold. Generally, there is no algorithmic way for such a manip- 
ulation. However, for polynomial functions all the related computations can 
be done algorithmically by means of Grobner bases, as we show in the next 
section. 
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3 Algorithm Description 



Here we describe an algorithm which, given a polynomial Lagrangian whose 
coefficients are rational numbers, computes the complete set of constraints 
and separates them into the first and second classes. The algorithm combines 
the above described Dirac method with the Grobner bases technique. By 
this reason we call it algorithm Dirac-Grobncr. All the below used concepts, 
definitions and constructive methods related to Grobner bases are explained. 



for instance, in textbooks E3,14 15 



At first we present the algorithm under assumption that a polynomial 
ideal generated by constraints is radical. This is true for most of real practical 
problems. Next, we indicate how to modify the algorithm to treat the most 
general (non-radical) case. 

Algorithm Dirac-Grobner 
Input: L{q,q), a polynomial Lagrangian {L e Q[q,q]) 

Output: <?i and $2, sets of the first and second class constraints, respectively. 

1. Computation of the canonical Hamiltonian and primary constraints: 

(a) Construct the polynomial set F — Uf^i{pi — dL/dqi} in variables 

p, q,q- 

(b) Compute the Grobner basis G of the ideal in ring Q [p, q, q] generated 
by F with respect to an ordering^ which eliminates q. Then compute 
the canonical Hamiltonian as the normal form of modulo G. 

(c) Find the set of primary constraint polynomials as G fl Q[p,q]. If 
<^ = 0, then stop since the system is regular. Otherwise, go to the 
next step. 

2. Computation of the complete set of constraints: 

(a) Take G — for the Grobner basis G of the ideal generated by ^ in 
Q[p,q] with respect to the ordering induced by that chosen at Step 
1(b). Fix this ordering in the sequel. 

(b) Construct the total Hamiltonian in form (|^) with multipliers Ua 
treated as symbolic constants (parameters). 

(c) For every element 4>a in <!> compute the normal form h of the Poisson 
bracket {iJt, ^q,} modulo G. li h ^ Q and no multipliers u/3 occur in 
it, then enlarge set <!> with h, and compute the Grobner basis G for 
the enlarged set. 

(d) If G = {1}, stop because the system is inconsistent. Otherwise, re- 
peat the previous step until the consistency condition (H) is satisfied 
for every element in ^ irrespective of multipliers Ua. This gives the 
complete set of constraints <1> — {0i, . . . , 0fc}. 



^ An elimination ordering which induced the degree-reverse-lexicographical one for 
monomials in p and q is heuristically best for efficiency reasons. 
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3. Separation of constraints into first and second classes: 

(a) Construct matrix M in by computing the normal forms of its 
elements modulo G, and determine rank r of M . If r = fc, stop with 
<l>^ = 0, = If J' = 0, stop with (Pi ^ <P and = 0- Otherwise, 
go to the next step. 

(b) Find a basis A = {ai, . . . , ak-r} of the null space (kernel) of the linear 
transformation defined by M. For every vector a in A construct a first 
class constraint as aa(l>a- Collect them in set 

(c) Construct {k — r) x k matrix {aj)a from components of vectors in A 
and find a basis B = {bi, . . . , br} of the null space of the correspond- 
ing linear transformation. For every vector b in B construct a second 
class constraint as ba(f>a- Collect them in set 'p2- 

The correctness of Steps 1, 2 and 3(a) of the algorithm is provided by the 
properties of Grobner bases [p^Ul^ , |l5| and by the following facts: (i) the def- 
inition (^) of the canonical Hamiltonian implies its independence of q on the 
primary constraint manifold (^) ; (ii) whenever a multiplier Ua in §) is differ- 
entiated when the Poisson bracket in is evaluated, the corresponding term 
vanishes on the constraint manifold. The correctness of Steps 3(b) and 3(c) 
follows from definition ( p^ of the first class constraints and the correctness 
of Step 3(a). The termination of algorithm Dirac-Grobner follows from the 
finiteness of the Grobner basis G which is constructed at Step 2(c). 

Now consider the most general case when the constraints obtained from 
(||) lead to a non-radical ideal. It should be noted that the ideal generated 
by the primary constraint polynomials (Step 1) is always radical. This is pro- 
vided by linearity of (^ in momenta. However, already the first secondary 
constraint added may destroy this property of the ideal. Therefore, the algo- 
rithm needs one more step, namely. Step 2(e), where the Grobner basis G of 
the radical ideal for the polynomial set <l> is computed. Next, every constraint 
polynomial in <P is replaced by its normal form modulo G. All the elements 
with zero normal forms are eliminated from the set. The extra step is also 
algorithmic. There are algorithms for construction of a basis, and, hence, a 
Grobner basis, of the radical of agiven ideal, which are built-in in some com- 
puter algebra systems (see [|l^,|lj,|5) for more details and references). One 
can also check the radical membership of h at Step 2(c) before its adding to 
(p. This check is easily done but in any case Step 2, for the correctness 

of Step 3, must end up with the radical sets ^ and G. 

We implemented algorithm Dirac-Grobner, as it presented above for the 
radical case, in Maple V Release 5. The implementation is relied on the built- 
in system facilities for computation and manipulation with Grobner bases and 
for linear algebra. Using our Maple code for different examples from physics 
and mechanics, we experimentally observed that in those infrequent cases 
when the constraint ideals are non-radical this can easily be detected from 
the structure of the output set. 
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4 Examples 



In this section we illustrate, by examples from physics and mechanics, the 
application of algorithm Dirac-Grobner. 

Example 1. SU{2) Yang- Mills mechanics in + 1 dimensional space-time 
This is a constrained physical model with gauge symmetry. The model La- 
grangian is given by L = = Xi + ge^jkUjXk (1 < ij, k < 

3). Here Xi and yi are the generalized coordinates and tensor eijk is anti- 
symmetric in its indices with £123 — 1. Respectively, the primary constraints 
and the canonical Hamiltonian are = and He ~ — ^ijkXjPkVi with the 
momenta given by p^ = dL/diji and pi = dL/dxi. The other constraints in 
the complete set computed by the algorithm are — CijkXjPk = 0, and all 
the six constraints found are of the first class. 

Example 2. Point particle of mass m moving on the surface of a sphere 
(rigid rotator). The movement is described by the regular Lagrangian io = 
^m^(g'i^ -I- -I- q3^)/2 = ^m'^q'^ with the externally imposed holonomic 
constraint = — 1 = 0. This system is equivalent to the singular La- 
grangian system L = Lq + A(/>, where A is an extra coordinate. There is the 
only primary constraint p\ — = dL/dX), and the canonical Hamilto- 
nian is He = ^m^p^ — X(p{q) {pi — dL/dqi). The complete set of constraint 
polynomials for the singular system contains four second class polynomials 
{pxi 4>{q),Piqi, 2mA Coming back to the initial regular system, the first 
and the last polynomials in the set must be omitted since they determine the 
extra dynamical variables. 

Example 3. Singular physical system with both first and second class con- 
straint^ The system Lagrangian is L = qi{q2 — qs) — qiq2- There are three 
primary constraint polynomials {pi + 92,^2 — 91,^3} ■ The canonical Hamil- 
tonian is He = qiqi- One more constraint polynomial qi is found by the 
Dirac-Grobner algorithm. The sets and ^2 of the first and second classes 
are {p2 + 91,^3} and {pi + 92 , 91 } , respectively. Note that this system has no 
physical degrees of freedom (c.f. 

Example 4- Inconsistent singular system L ~ ^q\+q2. There is the single 
primary constraint p2 = 0. The canonical Hamiltonian is He — p\/2 — q2. At 
Step 2(c) of algorithm Dirac-Grobner the inconsistency P2 — ^ occurs. The 
algorithm detects this inconsistency and stops. 

The above examples are rather small and can be treated by hand. With our 
Maple code we have already tried successfully much more nontrivial exam- 
ples. For instance, we computed and separated the constraints for the SU (2) 

* A.Burnel. Private communication. 
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Yang-Mills mechanics in 3 + 1 dimensional space-time |^ . Surprisingly, this 
computation took only a few seconds on an Pentium 100 personal computer 
though the model Lagrangian and the canonical Hamiltonian are rather cum- 
bersome polynomials of the 4th degree in 21 variables. 
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